Model evaluation: sensitivity and elasticity analysis for pelagic taxonomic units
1 Preamble
This is an Electronic Supplement to the manuscript Marques et al. “Quantifying Deepwater Horizon oil spill induced injury on pelagic cetaceans”.
The master file containing links to all supplementary files related to this paper is this Electronic Supplements master file.
If you make use of any of this material in your work, it would be appreciated if you would contact Tiago Marques to let him know.
2 Introduction
This working document describes how the injury results reported depend on the different population dynamics model parameters considered. All the parameters have uncertainty associated with them and it is only natural to ask what would change in the reported values for the injury metrics if the input parameters changed accordingly. We have 3 injury metrics:
- lost cetacean years (LCY),
- maximum proportional decrease (MPD), and
- years to recovery (YTR).
To do so we considered two types of sensitivity analysis:
Uncertainty analysis. The aim is to determine how much the uncertainty in each input parameter influences uncertainty in the resulting injury metrics. Each model parameter for which a distribution was available was allowed to vary by sampling from this distribution while keeping all the other parameters constant at their nominal value (i.e., the mean of their corresponding distributions). The relative contribution of uncertainty in each input parameter to overall uncertainty in injury metrics was thereby evaluated. This is akin to the sensitivity analysis undertaken by Schwacke et al. (2017 Supplemental Materials).
Elasticity analysis. The aim is to estimate the proportional change in the injury metrics caused by a proportional change in each model parameter. This allows us to quantify the “inherent sensitivity” of each model parameter. Elasticity is traditionally used in the population dynamics literature (Benton & Grant 199, Caswell 1981) to determine the proportional change in population growth rate to a proportional change in population vital rate; here we extend the concept to cover the injury metrics.
To contrast the two, imagine two model parameters for which uncertainty on inputs produces the same moderate level of uncertainty on the estimate of an injury metric. In one case this might be because there is considerable uncertainty on the input (i.e., the parameter is not known at all precisely) combined with low elasticity (i.e., little effect of this uncertainty on the injury metric); in the second case this might be because there is little uncertainty on the input combined with high elasticity. Undertaking both an uncertainty and an elasticity analysis allows us to distinguish between these scenarios.
The analyses reported in this document are based on files produced by sourcing the file RunAllPopSimsSensitivity. That source file produces for each species a set of result files, one for each parameter, that we then post-process here to obtain both uncertainty and elasticity values.
3 The model parameters
We provide a complete list of the quantities/parameters/statistics considered in population dynamics model, that might therefore have an impact in the injury metrics:
- Based on the models and results from Roberts et al., with distributions available via parametric bootstrap resampling; see also IniPopSizes for details.
- \(N_0\): initial population size
- \(p_e\): the proportion of the population exposed to oil
- Based on the Expert Elicitation exercise, see SI section ElectronicSupplement1_EE, for details
- \(\rho\): density dependence shape parameter; modeled as a shifted gamma distribution
- \(P_{recovery}\): proportion of exposed animals that recover to baseline survival and fecundity levels within their lifetime; modeled as a shifted beta distribution
- \(SR\): post-spill survival reduction; modeled as a shifted beta distribution
- Based on earlier published work on bottlenose dolphin data, then scaled for each taxonomic unit by Gestation Duration. For additional details regarding the bottlenose dolphin data, see also SI, Table S2.1 - Input parameters for population model in Schwacke et al. 2021 and Schwacke et al. 2017.
- \(S_m\) and \(S_f\): a vector of age and sex dependent survival for males and females, based on the Siler model parameters, see SurvivalReduction for details
- \(ASM\): age of sexual maturity; modeled as a gamma distribution; see AgeOfSexualMaturity for details.
- \(F_{max}\): maximum fecundity - modeled as a beta-PERT distribution
- \(F_{nom}\): nominal fecundity - modeled as a beta-PERT distribution
- \(Br_{Tt}\): baseline reproductive success - modeled as beta distribution (\(R_{baseline}\) in Schwacke et al. 2021)
- \(Por_{Tt}\): post-spill reproductive success - modeled as beta distribution (\(R_1\) in Schwacke et al. 2021)
All these quantities and their respective distributions/parameters are provided in file SpeciesDefinitionFile.xlsx.
A note regarding ASM. Given that this is a value for age in a model where age is incremented yearly, it needs to be an integer. As it turns out, the range of possible values of integers given the assumed distribution only covers, with non negligible probability, two to four different values depending on species. This is the reason why for the injury metric associated with ASM there are only a reduced number of different values, and that induces that the median and the lower/upper percentile of the injury distributions might sometimes coincide.
4 Sensitivity Analysis
# number of iterations used in the sensitivity analysis
runs4Sens <- 500
# pattern to look for the files relating to sensitivity analysis in the "InOutBySp/" folder
patt <- paste0(runs4Sens, "Sens")We want to estimate the sensitivity to the 11 parameters defined above. We have run 500 iterations for each parameter (and each species) to evaluate sensitivity, and then we derive the elasticity by post-processing these runs to evaluate how the injury metric changes by changing the parameter in the vicinity of its mean value. We address these in turn in the following two sub-sections.
4.1 Uncertainty
Consider a given parameter \(K\). Consider \(I\) as an example injury metric. The mean \(I\) from the 500 observations of a given injury metric is \(\bar{I}\). We represent the \(PC^k(I)\) as the proportional change in \(I\) for parameter \(K\).
We have represented in the plot below, for each injury metric and each parameter, the corresponding proportional change, in particular, the median, 2.5 and 97.5 percentiles of the distribution of proportional changes, across the 500 uncertainty iterations for each parameter \(K\)
\(PC^K_{j}(I)=\frac{I^K_{j}-\bar{I}}{\bar{I}},~j=1,2,...,500\)
where, for each parameter \(K\) for which uncertainty is being evaluated, \(K\) is allowed to change but all other parameters are kept at their nominal value.
4.2 Elasticity
In this analysis we used the values obtained under the uncertainty analysis - see previous sub-section - to model the percentage change in the 3 injury metrics (Lost Cetacean Years LCY, Years to Recovery YTR and Maximum Proportional Decrease MPD) as a function of a unit percentage change in the different input parameters. To do so, we used a generalized additive model (GAM), Gaussian response, and identity link. The reported results are virtually insensitive to the specific model chosen since the relations are all extremely smooth where it matters to evaluate elasticity, which is in the vicinity of the point estimate (here we considered the mean) of the parameter being evaluated. For a couple of injury metrics and parameters combinations the fit is far from ideal, in particular for YTR, but since this happens mostly when the elasticity is extremely small the effects are negligible. These problems originate because we are approximating change in a discrete injury metric (YTR) with a continuous function.
The corresponding results are shown below, where we report the percentage change on the 3 injury metrics when each parameter is varied by 1% around the mean (i.e. increased or decreased by 0.5%), while all the other parameters are kept at their mean value.
Here we calculate the injury metric percent change for each change in 1% on each of the 12 parameters. The next couple of code chunks do not output plots by default, but by removing the chunk options results='hide',fig.show='hide' the .Rmd will lead the code to show plots that provide insights to what is going on under the hood.
Note that for two of the inputs to the model are not scalars, namely
- \(S_{f,j}\) and \(S_{m,j}\), and
- \(P_{marked}\),
Since there is no univariate quantity to consider that the injury metric depends on the above procedure needs an extension, where we use a univariate measure to represent these quantities. We considered the weighted (weighted by the proportion in each age and sex class) average survival, which we refer to as meanSiler (we also look at these means by sex, namely meanSilerF for females and meanSilerM for males), and the average \(P_{marked}\), that we refer to as paPM, as the corresponding univariate quantities to evaluate the injury metrics against. Note that because the same mean value can be obtained with different combinations of the vectors, and different combinations of the corresponding vectors will have different impacts on the injury metric, this will lead to noise around the relationship we are trying to estimate.
5 Sensitivity results by taxonomic unit
The taxonomic units and the corresponding codes considered in this document are:
- Bwsp beaked whales Beaked whales spp
- Fatt pygmy killer whale Feresa attenuata
- Ggri Risso’s dolphin Grampus griseus
- Gmac short-finned pilot whale Globicephala macrorhynchus
- Kosp kogia Kogia sp.
- Pele melon-headed whale Peponocephala electra
- Pmac sperm whale Physeter macrocephalus
- Satt pantropical spotted dolphin Stenella attenuata
- Sbre rough-toothed dolphin Steno bredanensis
- Scly Clymene dolphin Stenella clymene
- Scoe striped dolphin Stenella coeruleoalba
- Sfro Atlantic spotted dolphin Stenella frontalis
- Slon spinner dolphin Stenella longirostris
- Ttro offshore bottlenose dolphins Tursiops truncatus
- Ttrs shelf bottlenose dolphins Tursiops truncatus
The results for each taxonomic unit are presented in turn below.
5.1 Beaked whales
5.1.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Beaked_whales"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.1.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 0.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Bwspsimres500SensascS.RData
## Elasticity analysis for Bwspsimres500SensASM.RData
## Elasticity analysis for Bwspsimres500SensBrTt.RData
## Elasticity analysis for Bwspsimres500SensFmax.RData
## Elasticity analysis for Bwspsimres500SensFnom.RData
## Elasticity analysis for Bwspsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 58 * EPS) is small (axis
## 2)
## Elasticity analysis for Bwspsimres500Senspe.RData
## Elasticity analysis for Bwspsimres500Sensper.RData
## Elasticity analysis for Bwspsimres500SensPorTt.RData
## Elasticity analysis for Bwspsimres500Sensrho.RData
## Elasticity analysis for Bwspsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.2 Pygmy killer whale
5.2.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Pygmy_killer_whale"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.2.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Fattsimres500SensascS.RData
## Elasticity analysis for Fattsimres500SensASM.RData
## Elasticity analysis for Fattsimres500SensBrTt.RData
## Elasticity analysis for Fattsimres500SensFmax.RData
## Elasticity analysis for Fattsimres500SensFnom.RData
## Elasticity analysis for Fattsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 56 * EPS) is small (axis
## 2)
## Elasticity analysis for Fattsimres500Senspe.RData
## Elasticity analysis for Fattsimres500Sensper.RData
## Elasticity analysis for Fattsimres500SensPorTt.RData
## Elasticity analysis for Fattsimres500Sensrho.RData
## Elasticity analysis for Fattsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.3 Risso’s dolphin
5.3.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Rissos_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.3.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Ggrisimres500SensascS.RData
## Elasticity analysis for Ggrisimres500SensASM.RData
## Elasticity analysis for Ggrisimres500SensBrTt.RData
## Elasticity analysis for Ggrisimres500SensFmax.RData
## Elasticity analysis for Ggrisimres500SensFnom.RData
## Elasticity analysis for Ggrisimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 48 * EPS) is small (axis
## 2)
## Elasticity analysis for Ggrisimres500Senspe.RData
## Elasticity analysis for Ggrisimres500Sensper.RData
## Elasticity analysis for Ggrisimres500SensPorTt.RData
## Elasticity analysis for Ggrisimres500Sensrho.RData
## Elasticity analysis for Ggrisimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.4 Short-finned pilot whale
5.4.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Pilot_whales"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.4.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Gmacsimres500SensascS.RData
## Elasticity analysis for Gmacsimres500SensASM.RData
## Elasticity analysis for Gmacsimres500SensBrTt.RData
## Elasticity analysis for Gmacsimres500SensFmax.RData
## Elasticity analysis for Gmacsimres500SensFnom.RData
## Elasticity analysis for Gmacsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 67 * EPS) is small (axis
## 2)
## Elasticity analysis for Gmacsimres500Senspe.RData
## Elasticity analysis for Gmacsimres500Sensper.RData
## Elasticity analysis for Gmacsimres500SensPorTt.RData
## Elasticity analysis for Gmacsimres500Sensrho.RData
## Elasticity analysis for Gmacsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.5 Kogia
5.5.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Kogia_whales"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.5.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Kospsimres500SensascS.RData
## Elasticity analysis for Kospsimres500SensASM.RData
## Elasticity analysis for Kospsimres500SensBrTt.RData
## Elasticity analysis for Kospsimres500SensFmax.RData
## Elasticity analysis for Kospsimres500SensFnom.RData
## Elasticity analysis for Kospsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 44 * EPS) is small (axis
## 2)
## Elasticity analysis for Kospsimres500Senspe.RData
## Elasticity analysis for Kospsimres500Sensper.RData
## Elasticity analysis for Kospsimres500SensPorTt.RData
## Elasticity analysis for Kospsimres500Sensrho.RData
## Elasticity analysis for Kospsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.6 Melon-headed whale
5.6.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Melon-headed_whale"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.6.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Pelesimres500SensascS.RData
## Elasticity analysis for Pelesimres500SensASM.RData
## Elasticity analysis for Pelesimres500SensBrTt.RData
## Elasticity analysis for Pelesimres500SensFmax.RData
## Elasticity analysis for Pelesimres500SensFnom.RData
## Elasticity analysis for Pelesimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 66 * EPS) is small (axis
## 2)
## Elasticity analysis for Pelesimres500Senspe.RData
## Elasticity analysis for Pelesimres500Sensper.RData
## Elasticity analysis for Pelesimres500SensPorTt.RData
## Elasticity analysis for Pelesimres500Sensrho.RData
## Elasticity analysis for Pelesimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.7 Sperm whale
5.7.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Sperm_whale"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.7.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Pmacsimres500SensascS.RData
## Elasticity analysis for Pmacsimres500SensASM.RData
## Elasticity analysis for Pmacsimres500SensBrTt.RData
## Elasticity analysis for Pmacsimres500SensFmax.RData
## Elasticity analysis for Pmacsimres500SensFnom.RData
## Elasticity analysis for Pmacsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 45 * EPS) is small (axis
## 2)
## Elasticity analysis for Pmacsimres500Senspe.RData
## Elasticity analysis for Pmacsimres500Sensper.RData
## Elasticity analysis for Pmacsimres500SensPorTt.RData
## Elasticity analysis for Pmacsimres500Sensrho.RData
## Elasticity analysis for Pmacsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
To further understand elasticity it might be useful to represent population trajectories under different extreme realizations of a given parameter. As an example, we show said trajectories for the mean Siler survival (i.e. mean Survival weighted for the proportion in each age class and sex, but not weighted for the proportion marked) below:
load(paste0("InOutBySp/",species,"/Pmacsimres500SensascS.RData"))
sr <- simres
par(mfrow=c(1,1),mar = c(4, 4, 2, 0.1))
nyears <- dim(sr)[2]
nsims <- dim(sr)[3]
ylims <- c(min(c(colSums(sr[, , 1:nsims, 2]), colSums(sr[, , 1:nsims, 1])))*0.9, max(c(colSums(sr[, , 1:nsims, 2]), colSums(sr[, , 1:nsims, 1])))*1.1)
plot(colSums(sr[, , 1, 2]), type = "n", ylim = ylims, xlab = "Years post spill", ylab = "Predicted population size", las = 1,main="Mean Siler")
iMAX <- which(pars2save$meanSiler==max(pars2save$meanSiler))[1]
iMIN <- which(pars2save$meanSiler==min(pars2save$meanSiler))[1]
# baseline under minimum value for parameter
lines(colSums(sr[, , iMIN, 2]), type = "l", lwd = 0.7, col = 2)
# baseline under maximum value for parameter
lines(colSums(sr[, , iMAX, 2]), type = "l", lwd = 1, col = 3)
# under oil minimum value for parameter
lines(colSums(sr[, , iMIN, 1]), type = "l", lwd = 0.7, col = 2)
# under oil maximum value for parameter
lines(colSums(sr[, , iMAX, 1]), type = "l", lwd = 1, col = 3)
# add a legend
legend("bottomright",inset=0.05,cex=1.2,col=2:3,lty=1,legend=c("min meanSiler","max meanSiler"))
In general terms, for each of the parameters, the relationship between what are the relative impacts in terms of uncertainty also reflect in similar patterns for elasticity. A noticeable difference is perhaps what happens for baseline survival, where the effect on the MPD presents an opposing pattern in elasticity. This pattern happens because while both YTR and LCY are naturally affected by the lower survival inducing longer recoveries, and hence much larger values for those injury metrics under low survival, the MPD is actually larger under high survival because the baseline scenario population size increases much more compared to post-spill scenario under the same conditions.
One might also want to identify why there is such variability in the MPD when considering the meanSiler. It is perhaps unanticipated why not so different “average” meanSiler survivals could lead to MPD ranging from 44% to 50%. Visualizing those below show us that in fact the trajectories are somewhat different. But if the meanSiler is similar, what causes the difference?
sr <- simres
par(mfrow=c(1,1),mar = c(4, 4, 2, 0.1))
nyears <- dim(sr)[2]
nsims <- dim(sr)[3]
ylims <- c(min(c(colSums(sr[, , 1:nsims, 2]), colSums(sr[, , 1:nsims, 1])))*0.9, max(c(colSums(sr[, , 1:nsims, 2]), colSums(sr[, , 1:nsims, 1])))*1.1)
plot(colSums(sr[, , 1, 2]), type = "n", ylim = ylims, xlab = "Years post spill", ylab = "Predicted population size", las = 1,main="Mean Siler")
iMAX <- which(injury$MPD==max(injury$MPD))[1]
iMIN <- which(injury$MPD==min(injury$MPD))[1]
# baseline under minimum value for parameter
lines(colSums(sr[, , iMIN, 2]), type = "l", lwd = 0.7, col = 2)
# baseline under maximum value for parameter
lines(colSums(sr[, , iMAX, 2]), type = "l", lwd = 1, col = 3)
# under oil minimum value for parameter
lines(colSums(sr[, , iMIN, 1]), type = "l", lwd = 0.7, col = 2)
# under oil maximum value for parameter
lines(colSums(sr[, , iMAX, 1]), type = "l", lwd = 1, col = 3)
# add a legend
legend("bottomright",inset=0.05,cex=1.2,col=2:3,lty=1,legend=c("min MPD","max MPD"))
Do we get a hint of what is going on if we plot these injury metrics against survival separated by sex?
par(mfrow = c(3, 3), mar = c(4, 4, 2.5, 0.5))
#for each possible x-axis value to evaluate elasticity against in the case of baseline survival
par.names<-c("meanSilerM","meanSilerF","meanSiler")
for(s in 1:3){
par.name <- par.names[s]
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=2) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par, k = 3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and addinmodel on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}
We can get a feeling for how female survival seems to be more relevant for the injury metrics if we investigate what induces the larger differences in MPD. Looking at the female and male survivals for the two iterations that lead to the largest differences in MPD, even though they have very similar meanSiler survivals of 0.939 (minimum MPD) and 0.91 (maximum MPD).
source("Functions/SilerFuns.R") # Siler model functions are here
naBND <- 61 # number of age classes for BND
agesBND <- 0:(naBND - 1)
#--------------------------------------------------------------------------
# Survival
#--------------------------------------------------------------------------
# reading in Siler parameters, for male and female BNDs
# These files were provided by LS via email on the Fri 5/22/2020 4:10 PM
pf <- read.csv("InputFiles/var_female 15May2020.csv")
pm <- read.csv("InputFiles/var_male 15May2020.csv")
# remove first column, which is just ID (not needed)
pf <- pf[, -1]
pm <- pm[, -1]
# get a sample from the Siler model for each sim - which is a posterior distribution
# number of Siler model draws
nS <- nrow(pf)
# objects to hold female and male realizations - Ttru
pxFs <- matrix(NA, nrow = length(agesBND), ncol = nS)
pxMs <- matrix(NA, nrow = length(agesBND), ncol = nS)
# for each observation of the posterior
for (s in 1:nS) {
# Survival Ttru
pxFs[, s] <- px(agesBND, pf[s, 1], pf[s, 2], pf[s, 3], pf[s, 4], pf[s, 5])
pxMs[, s] <- px(agesBND, pm[s, 1], pm[s, 2], pm[s, 3], pm[s, 4], pm[s, 5])
}
plot(agesBND,pxFs[,pars2save$iS[iMAX]],col="green",type="l",lty=1,ylab="Survival probability",xlab="Age")
points(agesBND,pxMs[,pars2save$iS[iMAX]],col="green",type="l",lty=2)
points(agesBND,pxFs[,pars2save$iS[iMIN]],col="red",type="l",lty=1)
points(agesBND,pxMs[,pars2save$iS[iMIN]],col="red",type="l",lty=2)
legend("bottomleft",inset=0.05,legend=c("Female max MPD","Male max MPD","Female min MPD","Male min MPD"),lty=c(1,2,1,2),col=c("green","green","red","red"))
#just taking a peak at these iteration parameters
#pars2save[c(iMIN,iMAX),]5.8 Pantropical spotted dolphin
5.8.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Pantropical_spotted_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.8.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Sattsimres500SensascS.RData
## Elasticity analysis for Sattsimres500SensASM.RData
## Elasticity analysis for Sattsimres500SensBrTt.RData
## Elasticity analysis for Sattsimres500SensFmax.RData
## Elasticity analysis for Sattsimres500SensFnom.RData
## Elasticity analysis for Sattsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 65 * EPS) is small (axis
## 2)
## Elasticity analysis for Sattsimres500Senspe.RData
## Elasticity analysis for Sattsimres500Sensper.RData
## Elasticity analysis for Sattsimres500SensPorTt.RData
## Elasticity analysis for Sattsimres500Sensrho.RData
## Elasticity analysis for Sattsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.9 Rough-toothed dolphin
5.9.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Rough-toothed_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.9.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Sbresimres500SensascS.RData
## Elasticity analysis for Sbresimres500SensASM.RData
## Elasticity analysis for Sbresimres500SensBrTt.RData
## Elasticity analysis for Sbresimres500SensFmax.RData
## Elasticity analysis for Sbresimres500SensFnom.RData
## Elasticity analysis for Sbresimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 68 * EPS) is small (axis
## 2)
## Elasticity analysis for Sbresimres500Senspe.RData
## Elasticity analysis for Sbresimres500Sensper.RData
## Elasticity analysis for Sbresimres500SensPorTt.RData
## Elasticity analysis for Sbresimres500Sensrho.RData
## Elasticity analysis for Sbresimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.10 Clymene dolphin
5.10.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Clymene_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.10.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Sclysimres500SensascS.RData
## Elasticity analysis for Sclysimres500SensASM.RData
## Elasticity analysis for Sclysimres500SensBrTt.RData
## Elasticity analysis for Sclysimres500SensFmax.RData
## Elasticity analysis for Sclysimres500SensFnom.RData
## Elasticity analysis for Sclysimres500SensN0.RData
## Elasticity analysis for Sclysimres500Senspe.RData
## Elasticity analysis for Sclysimres500Sensper.RData
## Elasticity analysis for Sclysimres500SensPorTt.RData
## Elasticity analysis for Sclysimres500Sensrho.RData
## Elasticity analysis for Sclysimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.11 Striped dolphin
5.11.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Striped_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.11.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Scoesimres500SensascS.RData
## Elasticity analysis for Scoesimres500SensASM.RData
## Elasticity analysis for Scoesimres500SensBrTt.RData
## Elasticity analysis for Scoesimres500SensFmax.RData
## Elasticity analysis for Scoesimres500SensFnom.RData
## Elasticity analysis for Scoesimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 46 * EPS) is small (axis
## 2)
## Elasticity analysis for Scoesimres500Senspe.RData
## Elasticity analysis for Scoesimres500Sensper.RData
## Elasticity analysis for Scoesimres500SensPorTt.RData
## Elasticity analysis for Scoesimres500Sensrho.RData
## Elasticity analysis for Scoesimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.12 Atlantic spotted dolphin
5.12.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Atlantic_spotted_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.12.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Sfrosimres500SensascS.RData
## Elasticity analysis for Sfrosimres500SensASM.RData
## Elasticity analysis for Sfrosimres500SensBrTt.RData
## Elasticity analysis for Sfrosimres500SensFmax.RData
## Elasticity analysis for Sfrosimres500SensFnom.RData
## Elasticity analysis for Sfrosimres500SensN0.RData
## Elasticity analysis for Sfrosimres500Senspe.RData
## Elasticity analysis for Sfrosimres500Sensper.RData
## Elasticity analysis for Sfrosimres500SensPorTt.RData
## Elasticity analysis for Sfrosimres500Sensrho.RData
## Elasticity analysis for Sfrosimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.13 Spinner dolphin
5.13.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Spinner_dolphin"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.13.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Slonsimres500SensascS.RData
## Elasticity analysis for Slonsimres500SensASM.RData
## Elasticity analysis for Slonsimres500SensBrTt.RData
## Elasticity analysis for Slonsimres500SensFmax.RData
## Elasticity analysis for Slonsimres500SensFnom.RData
## Elasticity analysis for Slonsimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 29 * EPS) is small (axis
## 2)
## Elasticity analysis for Slonsimres500Senspe.RData
## Elasticity analysis for Slonsimres500Sensper.RData
## Elasticity analysis for Slonsimres500SensPorTt.RData
## Elasticity analysis for Slonsimres500Sensrho.RData
## Elasticity analysis for Slonsimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.14 Offshore bottlenose dolphins
5.14.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Bottlenose_dolphin_oceanic"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.14.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Ttrosimres500SensascS.RData
## Elasticity analysis for Ttrosimres500SensASM.RData
## Elasticity analysis for Ttrosimres500SensBrTt.RData
## Elasticity analysis for Ttrosimres500SensFmax.RData
## Elasticity analysis for Ttrosimres500SensFnom.RData
## Elasticity analysis for Ttrosimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 40 * EPS) is small (axis
## 2)
## Elasticity analysis for Ttrosimres500Senspe.RData
## Elasticity analysis for Ttrosimres500Sensper.RData
## Elasticity analysis for Ttrosimres500SensPorTt.RData
## Elasticity analysis for Ttrosimres500Sensrho.RData
## Elasticity analysis for Ttrosimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
5.15 Shelf bottlenose dolphins
5.15.1 Uncertainty
We calculate the uncertainty measures per parameter and plot them:
species <- "Bottlenose_dolphin_shelf"
# get files names
files <- list.files(path = paste0("InOutBySp/",species), pattern = patt)
# check how many files
nfiles <- length(files)
# set up number of plots in figure
par(mfrow = c(1, 1), mar = c(5, 4, 2.5, 0.5))
plot(c(0.5, 11.5), c(-1, 1), type = "n", xlab = "", main = species, xaxt = "n", ylab = "Proportional change (95% CI)")
abline(h = 0, lty = 2)
for (i in 1:nfiles) {
# for each parameter
# load the sensitivity results
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by runs4Sens iterations data.frame
## a 3 column (one per injury metric) by nsims iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
# scale the injury by subtracting mean injury for each injury metric
sinjury <- scale(injury, scale = FALSE)
# divide by mean for each injury metric
sinjury[, 1] <- sinjury[, 1] / mean(injury[, 1], na.rm = TRUE)
sinjury[, 2] <- sinjury[, 2] / mean(injury[, 2], na.rm = TRUE)
sinjury[, 3] <- sinjury[, 3] / mean(injury[, 3], na.rm = TRUE)
# plot results
qsinjury <- apply(sinjury, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE)
segments(i - 0.2, qsinjury[1, 1], i - 0.2, qsinjury[3, 1], col = 4)
segments(i, qsinjury[1, 2], i, qsinjury[3, 2], col = 2)
segments(i + 0.2, qsinjury[1, 2], i + 0.2, qsinjury[3, 2], col = 1)
points(c(i - 0.2, i, i + 0.2), qsinjury[2, ], col = c(4, 2, 1), pch = c(16, 15, 18))
label <- par.name
if (par.name == "N0") label <- expression(N[0])
if (par.name == "br") label <- expression(R[baseline])
if (par.name == "rho") label <- expression(rho)
if (par.name == "per") label <- expression(P[recovery])
if (par.name == "Fmax") label <- expression(F[max])
if (par.name == "Fnom") label <- expression(F[nom])
if (par.name == "ascS") label <- expression(paste(S[f, j], " & ", S[m, j]))
if (par.name == "PM") label <- expression(P[marked])
axis(side = 1, at = i, labels = label, las = 2)
legend("topleft", inset = 0.05, legend = names(injury), col = c(4, 2, 1), lty = 1, pch = c(16, 15, 18))
}
5.15.2 Elasticity
We first define an object to hold the elasticity measure for each parameter.
metrics <- c("LCY", "YTR", "MPD")
nmetrics <- length(metrics)
elasticity <- data.frame(parID = 1:nfiles, parameter = NA, LCY = NA, YTR = NA, MPD = NA)Next we calculate the elasticity measures per parameter and plot them:
#```{r, results='hide', fig.show='hide'}
# set up number of plots in figure
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 0.5))
for (i in 1:nfiles) {
cat(paste0("Elasticity analysis for ", files[i]))
# this loads up a workspace that includes object injury
# a 3 column (one per injury metric) by number of iterations data.frame
load(paste0("InOutBySp/",species,"/", files[i]))
# for each parameter
# get name of parameter in files[i]
str <- files[i]
# parameter name start
startp <- str_locate(string = str, pattern = "Sens")[2] + 1
# parameters name end
endp <- str_locate(string = str, pattern = ".RD")[1] - 1
par.name <- substr(str, startp, endp)
elasticity$parameter[i] <- par.name
# select the parameter being evaluated - change parameter name later to avoid this tweak
if (par.name == "ascS") par.name <- "meanSiler"
if (par.name == "PM") par.name <- "paPM"
p2spar <- eval(parse(text = paste0("pars2save$", par.name)))
for (k in 1:nmetrics) {
# select the k_th injury (k = 1,2,3)
inj <- eval(parse(text = paste0("injury$", metrics[k])))
# plot results
plot(p2spar, inj, xlab = par.name, ylab = metrics[k])
dat4mod <- data.frame(par = p2spar, inj = injury[, k])
# if there are not enough unique values of the parameter being evaluated
# only happens if the parameter is ASM (2-3 unique values)
# force a linear model fit, else a relatively inflexible GAM
if (length(unique(inj))<=3) {
mod <- lm(inj ~ par, data = dat4mod)
} else {
mod <- gam(inj ~ s(par,k=3), data = dat4mod)
}
seqpars <- seq(min(p2spar), max(p2spar), length = 500)
preds <- predict(mod, newdata = data.frame(par = seqpars), type = "response")
# check model adequate for predictions by plotting relationship and add model on top
lines(seqpars, preds)
mpar <- mean(p2spar)
delta <- c(0.995, 1, 1.005)
nd <- data.frame(par = mpar * delta)
preds <- predict(mod, newdata = nd, type = "response")
m.metrics <- mean(inj)
abline(v = c(mpar * delta), col = c("red", "green", "red"))
abline(h = c(preds), col = c("red", "green", "red"))
# calculate elasticity measure for the m_th parameter and k_th injury metric
elas <- 100 * (preds[3] - preds[1]) / preds[2]
# save it in the right slot
elasticity[i, 3 + k] <- elas
# add number to top of the plot
mtext(text = round(elas,3), side = 3, cex = 0.8)
}
}## Elasticity analysis for Ttrssimres500SensascS.RData
## Elasticity analysis for Ttrssimres500SensASM.RData
## Elasticity analysis for Ttrssimres500SensBrTt.RData
## Elasticity analysis for Ttrssimres500SensFmax.RData
## Elasticity analysis for Ttrssimres500SensFnom.RData
## Elasticity analysis for Ttrssimres500SensN0.RData## Warning in plot.window(...): relative range of values ( 65 * EPS) is small (axis
## 2)
## Elasticity analysis for Ttrssimres500Senspe.RData
## Elasticity analysis for Ttrssimres500Sensper.RData
## Elasticity analysis for Ttrssimres500SensPorTt.RData
## Elasticity analysis for Ttrssimres500Sensrho.RData
## Elasticity analysis for Ttrssimres500SensSR.RData
We represent the elasticity measures per parameter here:
par(mfrow = c(1, 1), mar = c(4, 6, 2.5, 0.5))
barplot(t(as.matrix(elasticity[, 4:6])), col = c(4, 2, 1), beside = TRUE, ylab = "Percent change per 1% change in parameter", main = species)
for (i in 1:nrow(elasticity)) {
# set reasonable label
label <- switch(elasticity$parameter[i],
ASM = "ASM",
N0 = expression(N[0]),
br = expression(R[baseline]),
por = expression(R[1]),
rho = expression(rho),
per = expression(P[recovery]),
spos = expression(S[1]),
Fmax = expression(F[max]),
Fnom = expression(F[nom]),
ascS = expression(paste(S[f, j], " & ", S[m, j])),
PM = expression(P[marked]),
BrTt = "BrTt",
pe = expression(p[e]),
PorTt = "PorTt",
SR = "SR"
)
axis(side = 1, at = 2.5 + (i - 1) * 4, labels = label, las = 2)
}
legend("bottomleft", col = c(4, 2, 1), pt.bg = c(4, 2, 1), legend = metrics, inset = 0.2, cex = 1.2, pch = 22)
6 References
Benton, T. G. & Grant, A. 1999 Elasticity analysis as an important tool in evolutionary and population ecology Trends in Ecology & Evolution 14: 467-471
Caswell H. 2001 Matrix Population Models: Construction, Analysis, and Interpretation, 2nd ed. Sinauer Associates, Inc, Sunderland, Massachusetts.
Schwacke, L. H.; Thomas, L.; Wells, R. S.; McFee, W. E.; Hohn, A. A.; Mullin, K. D.; Zolman, E. S.; Quigley, B. M.; Rowles, T. K. & Schwacke, J. H. 2017. Quantifying injury to common bottlenose dolphins from the Deepwater Horizon oil spill using an age-, sex- and class-structured population model Endangered Species Research 33: 265-279